###################################################################################################
# Legislative Term Limits and Polarization ########################################################
# Michael Olson and Jon Rogowski ##################################################################
###################################################################################################

# This script generates results presented in Tables 1 and 2 in the text, and Tables A.1, B.1, B.2, 
# B.3, B.4, B.5, B.6, B.10, B.11, and Figures A.1, B.1 and B.2, in the supplementary materials. 

###################################################################################################
# working directory and packages ##################################################################
###################################################################################################

   #setwd("./termlimits_polarization_analysisdata")

#####################################################################################################
# Data ##############################################################################################
#####################################################################################################

  pol_dat <- read.csv("pol_dat16.csv")
  
  # make legislative professionalism logged for analyses
  
    pol_dat$legprofscore <- log(pol_dat$legprofscore)
    pol_dat$legprofscore_1993 <- log(pol_dat$legprofscore_1993)

#####################################################################################################
# Descriptive Statistics and Figures ################################################################
#####################################################################################################

# Figure A.1: Plot for State Term Limit Use #########################################################

  state <- sort(rep(state.abb,length(1993:2016)))
  year  <- rep(1993:2016,length(state.abb))
  plot_dat <- as.data.frame(cbind(state,year))
  plot_dat$year <- as.integer(as.character(plot_dat$year))
  plot_dat <- plot_dat[plot_dat$state!="NE",]
  
  plot_dat$house_term_limit_year <- rep(3000,nrow(plot_dat))

  plot_dat[plot_dat$st=="AZ","house_term_limit_year"] <- 2000
  plot_dat[plot_dat$st=="AR","house_term_limit_year"] <- 1998
  plot_dat[plot_dat$st=="CA","house_term_limit_year"] <- 1996
  plot_dat[plot_dat$st=="CO","house_term_limit_year"] <- 1998
  plot_dat[plot_dat$st=="FL","house_term_limit_year"] <- 2000
  plot_dat[plot_dat$st=="LA","house_term_limit_year"] <- 2007
  plot_dat[plot_dat$st=="ME","house_term_limit_year"] <- 1996
  plot_dat[plot_dat$st=="MI","house_term_limit_year"] <- 1998
  plot_dat[plot_dat$st=="MO","house_term_limit_year"] <- 2002
  plot_dat[plot_dat$st=="MT","house_term_limit_year"] <- 2000
  plot_dat[plot_dat$st=="NV","house_term_limit_year"] <- 2010
  plot_dat[plot_dat$st=="OH","house_term_limit_year"] <- 2000
  plot_dat[plot_dat$st=="OK","house_term_limit_year"] <- 2004
  plot_dat[plot_dat$st=="SD","house_term_limit_year"] <- 2000

  plot_dat$house_term_limit <- 0
  plot_dat[(plot_dat$year)>plot_dat$house_term_limit_year,"house_term_limit"] <- 1

  plot_dat1 <- aggregate(house_term_limit~year,data=plot_dat,FUN=mean)
  plot_dat2 <- plot_dat[plot_dat$year==1993 & plot_dat$house_term_limit_year<3000,c("state","house_term_limit_year")]
  colnames(plot_dat2) <- c("state","year")
  plot_dat2$year <- plot_dat2$year+1
  plot_dat2 <- merge(plot_dat2,plot_dat1,by=c("year"),all.x=TRUE)
  plot_dat2$state <- as.character(plot_dat2$state)

  plot_dat2 <- transform(plot_dat2, 
                         index = ave(state, year, 
                                     FUN = function(x) rank(x, ties.method = "first")))

  plot_dat21 <- plot_dat2[plot_dat2$index==1,]
  plot_dat21$index <- NULL ; colnames(plot_dat21)[2] <- "state1"
  plot_dat22 <- plot_dat2[plot_dat2$index==2,]
  plot_dat22$index <- NULL ; colnames(plot_dat22)[2] <- "state2"
  plot_dat23 <- plot_dat2[plot_dat2$index==3,]
  plot_dat23$index <- NULL ; colnames(plot_dat23)[2] <- "state3"
  plot_dat24 <- plot_dat2[plot_dat2$index==4,]
  plot_dat24$index <- NULL ; colnames(plot_dat24)[2] <- "state4"
  plot_dat25 <- plot_dat2[plot_dat2$index==5,]
  plot_dat25$index <- NULL ; colnames(plot_dat25)[2] <- "state5"

  plot_dat2 <- merge(plot_dat21,plot_dat22,by=c("year","house_term_limit"),all.x=TRUE)
  plot_dat2 <- merge(plot_dat2,plot_dat23,by=c("year","house_term_limit"),all.x=TRUE)
  plot_dat2 <- merge(plot_dat2,plot_dat24,by=c("year","house_term_limit"),all.x=TRUE)
  plot_dat2 <- merge(plot_dat2,plot_dat25,by=c("year","house_term_limit"),all.x=TRUE)
  
  plot_dat2[is.na(plot_dat2$state2),"state2"] <- ""
  plot_dat2[is.na(plot_dat2$state3),"state3"] <- ""
  plot_dat2[is.na(plot_dat2$state4),"state4"] <- ""
  plot_dat2[is.na(plot_dat2$state5),"state5"] <- ""
  
  plot_dat2$state_vector <- paste(plot_dat2$state1,plot_dat2$state2,plot_dat2$state3,plot_dat2$state4,plot_dat2$state5,sep=" ")
  plot_dat2$state_vector <- trimws(plot_dat2$state_vector)
  
  pdf(paste(output_path,"stepplot.pdf",sep=""),width=9,height=5)   
  step <- ggplot(data=plot_dat1,aes(x=year,y=house_term_limit))+
    geom_step(size=2)+
    geom_text(data=plot_dat2,aes(x=year,y=house_term_limit,label=state_vector),nudge_y=0.0075,vjust=0,size=4)+
    ylab("Share of State Lower Houses with Term Limits")+
    xlab("Year")+
    theme_minimal()
  print(step)
  dev.off()

# Table A.1: Summary Statistics of Key Variables ##################################################

  sumstat_sg <- stargazer(pol_dat[,c("l_diffs","house_term_limit","divided_gov",
                                     "legprofscore","abs_diff")],
                          summary.stat=c("mean","median","min","max","sd"),
                          covariate.labels=c("Legislative Polarization","Term Limits","Divided Gov."
                                             ,"ln(Leg. Professionalism)","Party Competitiveness"),
                          title="Summary Statistics of Key Variables",
                          table.layout ="-c-b-",table.placement = "!ht",
                          label="summary_stats")

  cat(sumstat_sg, sep = '\n', file = paste(output_path,"sumstat.tex",sep=""))


#####################################################################################################
# Analysis  #########################################################################################
#####################################################################################################
  
# Table 1: Fixed Effects OLS Estimates: State Legislative Polarization and Term Limits ##############

  pol_dat$term_limit_temp <- pol_dat$house_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$abs_diff

  pooled_nocov <- felm(l_diffs~term_limit_temp
                       |state+year|0|state
                       ,data=pol_dat)
  
  pooled_cov <- felm(l_diffs~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)
  
  pol_dat$temp_abs_diff <- pol_dat$hs_abs_diff
  
  house_nocov <- felm(h_diffs~term_limit_temp
                      |state+year|0|state
                      ,data=pol_dat)
  
  house_cov <- felm(h_diffs~term_limit_temp
                    +divided_gov
                    +legprofscore+temp_abs_diff
                    |state+year|0|state
                    ,data=pol_dat)

  pol_dat$term_limit_temp <- pol_dat$senate_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$sen_abs_diff
  
  senate_nocov <- felm(s_diffs~term_limit_temp
                       |state+year|0|state
                       ,data=pol_dat)
  
  senate_cov <- felm(s_diffs~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)

  panel_sg <- stargazer(pooled_nocov,pooled_cov,house_nocov,house_cov,senate_nocov,senate_cov,
                        add.lines = list(c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                         c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                         c("Projected $R^2$",round(summary(pooled_nocov)$P.r.squared,3),round(summary(pooled_cov)$P.r.squared,3)
                                           ,round(summary(house_nocov)$P.r.squared,3),round(summary(house_cov)$P.r.squared,3)
                                           ,round(summary(senate_nocov)$P.r.squared,3),round(summary(senate_cov)$P.r.squared,3))),
                        notes.append = FALSE,notes.label = "",
                        notes="\\parbox[t]{0.8\\textwidth}{\\scriptsize \\textit{Note}: Entries are linear regression coefficients with 
                        standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                        omit.stat = c("rsq", "f", "ser","adj.rsq"),font.size = "footnotesize",
                        star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                        dep.var.caption="Polarization",table.placement = "!ht",
                        dep.var.labels = c("Pooled","House","Senate"),
                        covariate.labels=c("Term Limits","Divided Gov.",
                                           "ln(Leg. Professionalism)","Party Competitiveness"),
                        label="baseline_panel_analysis",
                        table.layout ="-ld-#-t-as-n",
                        title="Term Limits and Polarization")
  
  cat(panel_sg, sep = '\n', file  = paste(output_path,"panel.tex",sep=""))


# Table 2: Asymmetric Polarization and Term Limits ################################################

  pol_dat$term_limit_temp <- pol_dat$house_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$abs_diff
  pol_dat$temp_dem        <- pol_dat$dem_leg_median
  pol_dat$temp_rep       <- pol_dat$rep_leg_median

  dem_pooled <- felm(temp_dem~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)
  
  rep_pooled <- felm(temp_rep~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)

  pol_dat$temp_abs_diff <- pol_dat$hs_abs_diff
  pol_dat$temp_dem        <- pol_dat$hou_dem
  pol_dat$temp_rep       <- pol_dat$hou_rep
  
  dem_house <- felm(temp_dem~term_limit_temp
                    +divided_gov
                    +legprofscore+temp_abs_diff
                    |state+year|0|state
                    ,data=pol_dat)
  
  rep_house <- felm(temp_rep~term_limit_temp
                    +divided_gov
                    +legprofscore+temp_abs_diff
                    |state+year|0|state
                    ,data=pol_dat)
  
  pol_dat$term_limit_temp <- pol_dat$senate_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$sen_abs_diff
  pol_dat$temp_dem        <- pol_dat$sen_dem
  pol_dat$temp_rep       <- pol_dat$sen_rep
  
  dem_senate <- felm(temp_dem~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)
  
  rep_senate <- felm(temp_rep~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)
  
  asym_sg <- stargazer(dem_pooled,dem_house,dem_senate,rep_pooled,rep_house,rep_senate,
                       add.lines = list(c("Chamber","Pooled","House","Senate","Pooled","House","Senate"),
                                        c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                        c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                        c("Projected $R^2$",round(summary(dem_pooled)$P.r.squared,3),round(summary(dem_house)$P.r.squared,3)
                                          ,round(summary(dem_senate)$P.r.squared,3),round(summary(rep_pooled)$P.r.squared,3)
                                          ,round(summary(rep_house)$P.r.squared,3),round(summary(rep_senate)$P.r.squared,3))),
                       notes.append = FALSE,notes.label = "",
                       notes="\\parbox[t]{0.8\\textwidth}{\\scriptsize \\textit{Note}: Entries are linear regression coefficients with 
                       standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                       omit.stat = c("rsq", "f", "ser","adj.rsq"),font.size = "footnotesize",
                       star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                       dep.var.labels=c("Democrats","Republicans"),
                       dep.var.caption = c("Party Medians"),table.placement = "!ht",
                       covariate.labels=c("Term Limits","Divided Gov.",
                                          "ln(Leg. Professionalism)","Party Competitiveness"),
                       label="asymmetric_polarization",
                       table.layout ="-ld-#-t-as-n",
                       title="Term Limits and Asymmetric Polarization")
  
  cat(asym_sg, sep = '\n', file = paste(output_path,"asym.tex",sep=""))

# Table B.1: Senate Polarization and Term Limits + House Term
  
  senate_plus_nocov <- felm(s_diffs~senate_tl_plus
                            |state+year|0|state
                            ,data=pol_dat)
  
  senate_plus_cov <- felm(s_diffs~senate_tl_plus
                          +divided_gov
                          +legprofscore+sen_abs_diff
                          |state+year|0|state
                          ,data=pol_dat)
  
  senate_panel_sg <- stargazer(senate_plus_nocov,senate_plus_cov,
                               add.lines = list(c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                                c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                                c("Projected $R^2$",round(summary(senate_plus_nocov)$P.r.squared,3),round(summary(senate_plus_cov)$P.r.squared,3))),
                               notes.append = FALSE,notes.label = "",
                               notes="\\parbox[t]{0.6\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with
                               standard errors clustered on states in parentheses. The main treatment variable is an indicator that takes a value of one
                               if senate term limits are in effect and at least one maximum house service period has passed since senate term limits were implemented.
                               $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                               omit.stat = c("rsq", "f", "ser","adj.rsq"),
                               star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                               dep.var.labels = c("Senate Polarization"),
                               covariate.labels=c("Senate Term Limits $+$ House Term","Divided Gov.",
                                                  "ln(Leg. Professionalism)","Party Competitiveness"),
                               label="senate_plus",table.placement = "!ht",
                               table.layout ="-ld-#-t-as-n",
                               title="Term Limits and Polarization: Senate Term Limits $+$ House Term")
  
  cat(senate_panel_sg, sep = '\n', file = paste(output_path,"senate_plus_panel.tex",sep=""))


# Table B.2: Robustness: State Legislative Polarization and Term-Limitedness ######################

  pol_dat$tled_temp <- pol_dat$term_limitedness
  pol_dat$term_limit_temp <- pol_dat$house_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$abs_diff
  
  
  pooled_tled <- felm(l_diffs~tled_temp
                      +divided_gov
                      +legprofscore+temp_abs_diff
                      |state+year|0|state
                      ,data=pol_dat)
  
  pooled_int_tled <- felm(l_diffs~tled_temp+term_limit_temp
                          +divided_gov
                          +legprofscore+temp_abs_diff
                          |state+year|0|state
                          ,data=pol_dat)
  
  pol_dat$tled_temp <- pol_dat$house_term_limitedness
  pol_dat$temp_abs_diff   <- pol_dat$hs_abs_diff
  
  
  house_tled <- felm(h_diffs~tled_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat)
  
  house_int_tled <- felm(h_diffs~tled_temp+term_limit_temp
                         +divided_gov
                         +legprofscore+temp_abs_diff
                         |state+year|0|state
                         ,data=pol_dat)
  
  pol_dat$tled_temp <- pol_dat$senate_term_limitedness
  pol_dat$term_limit_temp <- pol_dat$senate_term_limit  
  pol_dat$temp_abs_diff   <- pol_dat$sen_abs_diff
  
  
  senate_tled <- felm(s_diffs~tled_temp
                      +divided_gov
                      +legprofscore+temp_abs_diff
                      |state+year|0|state
                      ,data=pol_dat)
  
  senate_int_tled <- felm(s_diffs~tled_temp+term_limit_temp
                          +divided_gov
                          +legprofscore+temp_abs_diff
                          |state+year|0|state
                          ,data=pol_dat)
  
  tled_sg <- stargazer(pooled_tled,pooled_int_tled,house_tled,house_int_tled,senate_tled,senate_int_tled,
                       add.lines = list(c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                        c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                        c("Projected $R^2$",round(summary(pooled_tled)$P.r.squared,3),round(summary(pooled_int_tled)$P.r.squared,3)
                                          ,round(summary(house_tled)$P.r.squared,3),round(summary(house_int_tled)$P.r.squared,3)
                                          ,round(summary(senate_tled)$P.r.squared,3),round(summary(senate_int_tled)$P.r.squared,3))),
                       notes.append = FALSE,notes.label = "",
                       notes="\\parbox[t]{0.875\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                     standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                       omit.stat = c("rsq", "f", "ser","adj.rsq"),
                       star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,  
                       dep.var.caption = "Polarization",table.placement = "!ht",
                       dep.var.labels = c("Pooled","House","Senate"),
                       covariate.labels=c("Term Limitedness","Term Limits","Divided Gov."
                                          ,"ln(Leg. Professionalism)","Party Competitiveness"),
                       label="term_limitedness_robustness",
                       table.layout ="-ld-#-t-as-n",
                       title="Term Limits and Polarization: `Term Limitedness' Measure")
  
  cat(tled_sg, sep = '\n', file = paste(output_path,"tled.tex",sep=""))


# Table B.3: State Legislative Polarization and Term Limits: Additional Covariate Control #########

  pol_dat$term_limit_temp <- pol_dat$house_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$abs_diff
  
  ec_pooled_cov <- felm(l_diffs~term_limit_temp
                        +divided_gov
                        +legprofscore+temp_abs_diff
                        +gov_party
                        +logpop+pc_income
                        +unemploy
                        +pct_fb+gini
                        |state+year|0|state
                        ,data=pol_dat)
  
  pol_dat$temp_abs_diff <- pol_dat$hs_abs_diff
  
  ec_house_cov <- felm(h_diffs~term_limit_temp
                       +divided_gov
                       +legprofscore+temp_abs_diff
                       +gov_party
                       +logpop+pc_income
                       +unemploy
                       +pct_fb+gini
                       |state+year|0|state
                       ,data=pol_dat)
  
  pol_dat$term_limit_temp <- pol_dat$senate_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$sen_abs_diff
  
  ec_senate_cov <- felm(s_diffs~term_limit_temp
                        +divided_gov
                        +legprofscore+temp_abs_diff
                        +gov_party
                        +logpop+pc_income
                        +unemploy
                        +pct_fb+gini
                        |state+year|0|state
                        ,data=pol_dat)
  
  extra_controls_sg <- stargazer(ec_pooled_cov,ec_house_cov,ec_senate_cov,
                                 add.lines = list(c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                                  c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                                  c("Projected $R^2$",round(summary(ec_pooled_cov)$P.r.squared,3)
                                                    ,round(summary(ec_house_cov)$P.r.squared,3)
                                                    ,round(summary(ec_senate_cov)$P.r.squared,3))),
                                 notes.append = FALSE,notes.label = "",
                                 notes="\\parbox[t]{0.525\\textwidth}{\\scriptsize \\textit{Note}: Entries are linear regression coefficients with 
                                 standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                                 omit.stat = c("rsq", "f", "ser","adj.rsq"),font.size = "footnotesize",
                                 star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                                 dep.var.caption="Polarization",table.placement = "!ht",
                                 dep.var.labels = c("Pooled","House","Senate"),
                                 covariate.labels=c("Term Limits","Divided Gov.","ln(Leg. Professionalism)","Party Competitiveness","Democratic Governor"
                                                    ,"ln(Population)","Per Capita Income","Unemployment Rate"
                                                    ,"Percent Foreign Born","State Gini Coefficient"),
                                 label="extra_covariates",
                                 table.layout ="-ld-#-t-as-n",
                                 title="Term Limits and Polarization: Additional Covariate Control")
  
  cat(extra_controls_sg, sep = '\n', file = paste(output_path,"extra_controls.tex",sep=""))


# Table B.4: State Legislative Polarization and Term Limits: Pre-Treatment Covariate Control ######
  
  ag_base_cov1 <- lm(l_diffs~house_term_limit
                     +(divided_gov_1993+legprofscore_1993+abs_diff_1993+gov_party_1993)*factor(year)
                     +factor(state)+factor(year)-1
                     ,data=pol_dat)
  ag_base_cov1_vcov <- cluster.vcov(ag_base_cov1,pol_dat[,"state"])
  
  ag_base_cov2 <- lm(l_diffs~house_term_limit
                     +(divided_gov_1993+legprofscore_1993+abs_diff_1993
                       +gov_party_1993+logpop_1993+pc_income_1993
                       +unemploy_1993)*factor(year)
                     +factor(state)+factor(year)-1
                     ,data=pol_dat)
  ag_base_cov2_vcov <- cluster.vcov(ag_base_cov2,pol_dat[,"state"])
  
  ag_base_cov3 <- lm(l_diffs~house_term_limit
                     +(divided_gov_1993+gov_party_1993+legprofscore_1993+abs_diff_1993
                       +logpop_1993+pc_income_1993
                       +unemploy_1993
                       +pct_fb_1993+gini_1993)*factor(year)
                     +factor(state)+factor(year)-1
                     ,data=pol_dat)
  ag_base_cov3_vcov <- cluster.vcov(ag_base_cov3,pol_dat[,"state"])
  
  ptcovs_sg <- stargazer(ag_base_cov1,ag_base_cov2,ag_base_cov3,
                         se=list(sqrt(diag(ag_base_cov1_vcov)),
                                 sqrt(diag(ag_base_cov2_vcov)),sqrt(diag(ag_base_cov3_vcov))),
                         keep=c("house_term_limit"),
                         add.lines = list(c("Covariates","Institutional","$+$ Demographic","$+$ MPR"),
                                          c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark"),
                                          c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark")),
                         notes.append = FALSE,notes.label = "",
                         notes="\\parbox[t]{0.7\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                         standard errors clustered on states in parentheses. Models control for covariates held at 1993 level and interacted with year indicators.
                         'Institutional' model includes controls for the base rates of divided government,
                         logged legislative professionalism, party competitiveness, and the party of the governor. 'Demographic' model adds to these
                         controls for logged state population, per-capita income, and the unemployment rate. 'MPR' model adds to these controls for
                         percent of population that is foreign-born and the state's gini coefficient.
                         $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                         omit.stat = c("rsq", "f", "ser","adj.rsq"),
                         star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                         dep.var.labels = c("Legislative Polarization"),
                         covariate.labels=c("Term Limits","Divided Gov.","Democratic Governor","ln(Leg. Professionalism)","Democratic Seat Share","Party Competitiveness"
                                            ,"ln(Population)","Per Capita Income","GOP Presidential State Share","Unemployment Rate"
                                            ,"Percent Foreign Born","State Gini Coefficient"),
                         label="baseline_robustness",
                         table.layout ="-ld-#-t-as-n",table.placement = "!ht",
                         title="Term Limits and Polarization: Pre-Treatment Covariate Control")
  
  cat(ptcovs_sg, sep = '\n', file = paste(output_path,"ptcovs.tex",sep=""))


# Table B.5: Robustness: Lagged Dependent Variable Model ############################################################################

  pol_dat$term_limit_temp <- pol_dat$house_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$abs_diff
  pol_dat$temp_ldv <- pol_dat$l_diffs_lag
  
  pooled_ldv_nocov <- felm(l_diffs~term_limit_temp
                           +temp_ldv
                           |year|0|state
                           ,data=pol_dat)
  
  pooled_ldv_cov <- felm(l_diffs~term_limit_temp
                         +divided_gov
                         +legprofscore+temp_abs_diff
                         +temp_ldv
                         |year|0|state
                         ,data=pol_dat)
  
  pol_dat$temp_abs_diff <- pol_dat$hs_abs_diff
  pol_dat$temp_ldv <- pol_dat$h_diffs_lag
  
  house_ldv_nocov <- felm(h_diffs~term_limit_temp
                          +temp_ldv
                          |year|0|state
                          ,data=pol_dat)
  
  house_ldv_cov <- felm(h_diffs~term_limit_temp
                        +divided_gov
                        +legprofscore+temp_abs_diff
                        +temp_ldv
                        |year|0|state
                        ,data=pol_dat)
  
  pol_dat$term_limit_temp <- pol_dat$senate_term_limit
  pol_dat$temp_abs_diff   <- pol_dat$sen_abs_diff
  pol_dat$temp_ldv <- pol_dat$s_diffs_lag
  
  senate_ldv_nocov <- felm(s_diffs~term_limit_temp
                           +temp_ldv
                           |year|0|state
                           ,data=pol_dat)
  
  senate_ldv_cov <- felm(s_diffs~term_limit_temp
                         +divided_gov
                         +legprofscore+temp_abs_diff
                         +temp_ldv
                         |year|0|state
                         ,data=pol_dat)
   
  ldv_sg <- stargazer(pooled_ldv_nocov,pooled_ldv_cov,
                      house_ldv_nocov,house_ldv_cov,
                      senate_ldv_nocov,senate_ldv_cov,
                      keep=c("term_limit","divided_gov","gov_party","legprofscore","dem_leg_share"
                             ,"temp_abs_diff","temp_ldv"),
                      add.lines = list(c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                       c("Projected $R^2$",round(summary(pooled_ldv_nocov)$P.r.squared,3),round(summary(pooled_ldv_cov)$P.r.squared,3),
                                         round(summary(house_ldv_nocov)$P.r.squared,3),round(summary(house_ldv_cov)$P.r.squared,3),
                                         round(summary(senate_ldv_nocov)$P.r.squared,3),round(summary(senate_ldv_cov)$P.r.squared,3))),
                      notes.append = FALSE,notes.label = "",
                      star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                      notes="\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with
                      standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                      omit.stat = c("rsq", "f", "ser","adj.rsq"),
                      dep.var.caption = c("Legislative Polarization"),
                      dep.var.labels = c("Pooled","House","Senate"),
                      covariate.labels=c("Term Limits","Divided Gov.",
                                         "ln(Leg. Professionalism)",
                                         "Party Competitiveness","Lagged Polarization"),
                      label="robustness_lag",
                      table.layout ="-ld-#-t-as-n",table.placement = "!ht",
                      title="Term Limits and Polarization: Lagged Dependent Variable Model")
  
  cat(ldv_sg, sep = '\n', file = paste(output_path,"ldv.tex",sep=""))
 

# Table B.6: Matches Sample Fixed Effects Model: State Legislative Polarization and Term Limits ###

  match_dat <- read.csv("pol_dat_formatch16.csv")
  match_dat$legprofscore <- log(match_dat$legprofscore)
  match_dat <- dplyr::arrange(match_dat,state,year)
  
  match_dat <- match_dat[match_dat$state!="NE",] 

  dat_list <- list()
  
  for(i in 1993:2016){
    
    dat <- match_dat[match_dat$year==i,c("divided_gov","legprofscore","abs_diff")]
    colnames(dat) <- paste(colnames(dat),i)
    
    dat_list[[i-1992]] <- dat
    
  }
  
  pol_dat_match <- do.call("bind_cols",dat_list)
  
  pol_dat_match <- as.data.frame(cbind(state.abb[state.abb!="NE"],pol_dat_match))
  colnames(pol_dat_match)[1] <- c("state")
  
  pol_dat_match$treated <- 0
  pol_dat_match[is.element(pol_dat_match$state,unique(pol_dat[pol_dat$house_term_limit==1,"st"])),"treated"] <- 1
  
  match_out <- Match(Tr=pol_dat_match$treated,
                     X=pol_dat_match[,2:which(colnames(pol_dat_match)=="abs_diff 1996")],
                     Weight=1,M=1)
 
  match_states <- pol_dat_match$state[c(match_out$index.control,match_out$index.treated)]
  
  pd_matched <- pol_dat[is.element(pol_dat$state,match_states),]

  match_nocov_11 <- felm(l_diffs~house_term_limit
                         |state+year|0|state
                         ,data=pd_matched)

  match_cov3_11 <- felm(l_diffs~house_term_limit
                        +divided_gov
                        +legprofscore+abs_diff
                        |state+year|0|state
                        ,data=pd_matched)


  match_out2 <- Match(Tr=pol_dat_match$treated,
                      X=pol_dat_match[,2:which(colnames(pol_dat_match)=="abs_diff 1996")],
                      Weight=1,M=2)

  match_states2 <- pol_dat_match$state[c(match_out2$index.control,match_out2$index.treated)]
  
  pd_matched2 <- pol_dat[is.element(pol_dat$state,match_states2),]
  
  match_nocov_21 <- felm(l_diffs~house_term_limit
                         |state+year|0|state
                         ,data=pd_matched2)
  
  match_cov3_21 <- felm(l_diffs~house_term_limit
                        +divided_gov
                        +legprofscore+abs_diff
                        |state+year|0|state
                        ,data=pd_matched2)
  
  match_sg <- stargazer(match_nocov_11, match_cov3_11, match_nocov_21, match_cov3_21,
                        add.lines = list(c("Matching","One-to-One","One-to-One","Two-to-One","Two-to-One"),
                                         c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                         c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                         c("Projected $R^2$",round(summary(match_nocov_11)$P.r.squared,3),round(summary(match_cov3_11)$P.r.squared,3)
                                           ,round(summary(match_nocov_21)$P.r.squared,3),round(summary(match_cov3_21)$P.r.squared,3))),
                        notes.append = FALSE,notes.label = "",
                        notes="\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                        standard errors clustered on states in parentheses. Sample is matched using inverse-variance weighting nearest-neighbor matching
                        on the basis of pre-treatment (1993-1996) covariates 
                        using the \\texttt{Matching} package in \\texttt{R}.
                        $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                        omit.stat = c("rsq", "f", "ser","adj.rsq"),
                        star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                        dep.var.labels = c("Legislative Polarization"),table.placement = "!ht",
                        covariate.labels=c("Term Limits","Divided Gov.",
                                           "ln(Leg. Professionalism)","Party Competitiveness"),
                        label="matched_panel_analysis",
                        table.layout ="-ld-#-t-as-n",
                        title="Term Limits and Polarization: Matched Sample")
  
  cat(match_sg, sep = '\n', file = paste(output_path,"match.tex",sep=""))

# Figure B.1: Coefficients when Successively Dropping States ######################################
  
  coef <- rep(NA,length(unique(pol_dat$state))); se <- coef; dropped <- coef
  pol_dat$state <- as.character(pol_dat$state)
  
  for(i in 1:length(unique(pol_dat$state))){
    
    loop <- pol_dat[pol_dat$state!=unique(pol_dat$state)[i],]
    
    loopreg <- felm(l_diffs~house_term_limit
                    +divided_gov
                    +legprofscore+abs_diff
                    |state+year|0|state
                    ,data=loop)
    
    coef[i] <- coef(loopreg)["house_term_limit"]
    se[i]   <- sqrt(loopreg$clustervcv["house_term_limit","house_term_limit"])
    dropped[i] <- unique(pol_dat$state)[i]
    
  }
  
  min(coef); max(coef)
  
  dropstate <- as.data.frame(cbind(dropped,coef,se))
  dropstate$coef <- as.numeric(as.character(dropstate$coef))
  dropstate$se <- as.numeric(as.character(dropstate$se))
  dropstate$upper95 <- dropstate$coef+2*dropstate$se; dropstate$lower95 <- dropstate$coef-2*dropstate$se
  dropstate$upper90 <- dropstate$coef+1.65*dropstate$se; dropstate$lower95 <- dropstate$coef-1.65*dropstate$se
  
  pdf(paste(output_path,"dropstate.pdf",sep=""),width=10,height=4)
  dropstate <- ggplot(data=dropstate,aes(x=dropped,y=coef,ymax=upper95,ymin=lower95))+
    geom_point()+
    geom_errorbar()+
    geom_hline(yintercept=0,colour="indianred",linetype=2,size=1.5)+
    ylab("Coefficient on 'Term Limits'")+
    xlab("Dropped State")+
    theme_minimal()
  print(dropstate)
  dev.off()

# Figure B.2: Coefficients when Successively Dropping Years #######################################
  
  coef2 <- rep(NA,length(unique(pol_dat$year))); se2 <- coef2; dropped2 <- coef2
  
  for(i in 1:length(unique(pol_dat$year))){
    
    loop2 <- pol_dat[pol_dat$year!=unique(pol_dat$year)[i],]
    
    loopreg2 <- felm(l_diffs~house_term_limit
                     +divided_gov
                     +legprofscore+abs_diff
                     |state+year|0|state
                     ,data=loop2)
    
    coef2[i] <- coef(loopreg2)["house_term_limit"]
    se2[i]   <- sqrt(loopreg2$clustervcv["house_term_limit","house_term_limit"])
    dropped2[i] <- unique(pol_dat$year)[i]
    
    print(i) 
    
  }
  
  min(coef2); max(coef2)
  
  dropyear <- as.data.frame(cbind(dropped2,coef2,se2))
  colnames(dropyear) <- c("dropped","coef","se")
  dropyear$coef <- as.numeric(as.character(dropyear$coef))
  dropyear$se <- as.numeric(as.character(dropyear$se))
  dropyear$upper95 <- dropyear$coef+2*dropyear$se; dropyear$lower95 <- dropyear$coef-2*dropyear$se
  dropyear$upper90 <- dropyear$coef+1.65*dropyear$se; dropyear$lower95 <- dropyear$coef-1.65*dropyear$se
  
  pdf(paste(output_path,"dropyear.pdf",sep=""),width=10,height=4)
  dropyear <- ggplot(data=dropyear,aes(x=dropped,y=coef,ymax=upper95,ymin=lower95))+
    geom_point()+
    geom_errorbar()+
    geom_hline(yintercept=0,colour="indianred",linetype=2,size=1.5)+
    ylab("Coefficient on 'Term Limits'")+
    xlab("Dropped Year")+
    theme_minimal()
  print(dropyear)
  dev.off()
  
# Table B.10: Fixed Effects OLS Estimates, with Nebraska ##############
  
  pol_dat_ne <- read.csv("pol_dat_wNE_16.csv")
  
  pol_dat_ne$term_limit_temp <- pol_dat_ne$house_term_limit
  pol_dat_ne$temp_abs_diff   <- pol_dat_ne$abs_diff
  
  pooled_ne_nocov <- felm(l_diffs~term_limit_temp
                       |state+year|0|state
                       ,data=pol_dat_ne)
  
  pooled_ne_cov <- felm(l_diffs~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat_ne)
  
  pol_dat_ne$term_limit_temp <- pol_dat_ne$senate_term_limit
  pol_dat_ne$temp_abs_diff   <- pol_dat_ne$sen_abs_diff
  
  senate_ne_nocov <- felm(s_diffs~term_limit_temp
                       |state+year|0|state
                       ,data=pol_dat_ne)
  
  senate_ne_cov <- felm(s_diffs~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat_ne)
  
  panel_ne_sg <- stargazer(pooled_ne_nocov,pooled_ne_cov,
                           senate_ne_nocov,senate_ne_cov,
                        add.lines = list(c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                         c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                         c("Projected $R^2$",round(summary(pooled_ne_nocov)$P.r.squared,3),round(summary(pooled_ne_cov)$P.r.squared,3)
                                           ,round(summary(senate_ne_nocov)$P.r.squared,3),round(summary(senate_ne_cov)$P.r.squared,3))),
                        notes.append = FALSE,notes.label = "",
                        notes="\\parbox[t]{0.675\\textwidth}{\\scriptsize \\textit{Note}: Entries are linear regression coefficients with 
                        standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                        omit.stat = c("rsq", "f", "ser","adj.rsq"),
                        star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                        dep.var.caption="Polarization",table.placement = "!ht",
                        dep.var.labels = c("Pooled","Senate"),
                        covariate.labels=c("Term Limits","Divided Gov.",
                                           "ln(Leg. Professionalism)","Party Competitiveness"),
                        label="baseline_panel_analysis_ne",
                        table.layout ="-ld-#-t-as-n",
                        title="Term Limits and Polarization: Including Nebraska")
  
  cat(panel_ne_sg, sep = '\n', file  = paste(output_path,"panel_ne.tex",sep=""))
  
  
# Table B.11: Asymmetric Polarization and Term Limits, with Nebraska ################################################
  
  pol_dat_ne$term_limit_temp <- pol_dat_ne$house_term_limit
  pol_dat_ne$temp_abs_diff   <- pol_dat_ne$abs_diff
  pol_dat_ne$temp_dem        <- pol_dat_ne$dem_leg_median
  pol_dat_ne$temp_rep       <- pol_dat_ne$rep_leg_median
  
  dem_ne_pooled <- felm(temp_dem~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat_ne)
  
  rep_ne_pooled <- felm(temp_rep~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat_ne)
  
  pol_dat_ne$term_limit_temp <- pol_dat_ne$senate_term_limit
  pol_dat_ne$temp_abs_diff   <- pol_dat_ne$sen_abs_diff
  pol_dat_ne$temp_dem        <- pol_dat_ne$sen_dem
  pol_dat_ne$temp_rep       <- pol_dat_ne$sen_rep
  
  dem_ne_senate <- felm(temp_dem~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat_ne)
  
  rep_ne_senate <- felm(temp_rep~term_limit_temp
                     +divided_gov
                     +legprofscore+temp_abs_diff
                     |state+year|0|state
                     ,data=pol_dat_ne)
  
  asym_ne_sg <- stargazer(dem_ne_pooled,
                       dem_ne_senate,rep_ne_pooled,
                       rep_ne_senate,
                       add.lines = list(c("Chamber","Pooled","Senate","Pooled","Senate"),
                                        c("State Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                        c("Year Fixed Effects","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                                        c("Projected $R^2$",round(summary(dem_ne_pooled)$P.r.squared,3)
                                          ,round(summary(dem_ne_senate)$P.r.squared,3),round(summary(rep_ne_pooled)$P.r.squared,3)
                                          ,round(summary(rep_ne_senate)$P.r.squared,3))),
                       notes.append = FALSE,notes.label = "",
                       notes="\\parbox[t]{0.675\\textwidth}{\\scriptsize \\textit{Note}: Entries are linear regression coefficients with 
                       standard errors clustered on states in parentheses.  $^{**}$p$<$0.05, $^*$p$<$0.10 (two-tailed test).}",
                       omit.stat = c("rsq", "f", "ser","adj.rsq"),
                       star.char=c("*","**"),star.cutoffs = c(0.10,0.05),digits=3,digits.extra=0,
                       dep.var.labels=c("Democrats","Republicans"),
                       dep.var.caption = c("Party Medians"),table.placement = "!ht",
                       covariate.labels=c("Term Limits","Divided Gov.",
                                          "ln(Leg. Professionalism)","Party Competitiveness"),
                       label="asymmetric_polarization_ne",
                       table.layout ="-ld-#-t-as-n",
                       title="Term Limits and Asymmetric Polarization: Including Nebraska")
  
  cat(asym_ne_sg, sep = '\n', file = paste(output_path,"asym_ne.tex",sep=""))


# within-state SD (referenced in Paragraph 1 of "Panel Evidence" section, and Footnote 15)

  sdag <- aggregate(l_diffs~state,FUN=sd,data=pol_dat)
  
  mean(sdag$l_diffs); median(sdag$l_diffs)
  
  tl_states <- unique(pol_dat$state[pol_dat$house_term_limit==1])
  
  sdag <- aggregate(l_diffs~state,FUN=sd,data=pol_dat[!is.element(pol_dat$state,tl_states),])
  
  mean(sdag$l_diffs); median(sdag$l_diffs)

# reset working directory
  
  setwd(storewd)
  